#! /usr/bin/env python

import math

def poly1(Xo, Xn, x):

    # 21st-order Polynomial

    P = [ -0.088098528242219,
          0.000000000115645,
          0.973720575237457,
          -0.000000001300817,
          -4.897241716183775,
          0.000000006689319,
          14.800552740472247,
          -0.000000020812281,
          -29.885731490355798,
          0.000000043721985,
          42.383401012952483,
          -0.000000065654557,
          -43.168278794320933,
          0.000000073100692,
          31.715470115427561,
          -0.000000062681098,
          -16.650621789606145,
          0.000000042328202,
          6.166896940467320,
          -0.000000021784738,
          -1.850069065848198,
          0.500000006277647]

    dimP = len(P)
    medx = 0.5*(Xo + Xn)
    a0 = 0.2*(Xn - Xo)
    r = abs(medx - x)/(0.5 * a0) - 1.0

    # Initial Pulse
    y = 0.0
    if r <= 1.0:
        for n in range(dimP):
            y += P[n] * math.pow(r, dimP - (n + 1))

    print("{} {}".format(x, y))

def poly2(x):

    Xo = -0.5
    Xn = 0.5
    medx = 0.5*(Xo + Xn)
    a0 = 0.2*(Xn - Xo)
    r = abs(medx - x)/(0.5 * a0) - 1.0

    c21 = -0.088098528242219
    c20 =  0.000000000115645
    c19 =  0.973720575237457
    c18 = -0.000000001300817
    c17 = -4.897241716183775
    c16 =  0.000000006689319
    c15 =  14.800552740472247
    c14 = -0.000000020812281
    c13 = -29.885731490355798
    c12 =  0.000000043721985
    c11 =  42.383401012952483
    c10 = -0.000000065654557
    c9  = -43.168278794320933
    c8  =  0.000000073100692
    c7  =  31.715470115427561
    c6  = -0.000000062681098
    c5  = -16.650621789606145
    c4  =  0.000000042328202
    c3  =  6.166896940467320
    c2  = -0.000000021784738
    c1  = -1.850069065848198
    c0  =  0.500000006277647

    value = 0.0
    if r <= 1:
        value = c0  * pow(r, 0) + \
                c1  * pow(r, 1) + \
                c2  * pow(r, 2) + \
                c3  * pow(r, 3) + \
                c4  * pow(r, 4) + \
                c5  * pow(r, 5) + \
                c6  * pow(r, 6) + \
                c7  * pow(r, 7) + \
                c8  * pow(r, 8) + \
                c9  * pow(r, 9) + \
                c10 * pow(r, 10) + \
                c11 * pow(r, 11) + \
                c12 * pow(r, 12) + \
                c13 * pow(r, 13) + \
                c14 * pow(r, 14) + \
                c15 * pow(r, 15) + \
                c16 * pow(r, 16) + \
                c17 * pow(r, 17) + \
                c18 * pow(r, 18) + \
                c19 * pow(r, 19) + \
                c20 * pow(r, 20) + \
                c21 * pow(r, 21)

    print("{} {}".format(x, value))

if __name__ == "__main__":

    x0 = -0.5
    xn = 0.5
    n_pts = 1000
    vals = [x*(xn - x0)/1000 + x0 for x in range(n_pts+1)]

    for val in vals:
        poly2(val)

